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Abstract 

An algorithm is given in this paper for the computation of dynamically equivalent weakly 
reversible realizations with the maximal number of reactions, for chemical reaction networks 
(CRNs) with mass action kinetics. The original problem statement can be traced back at least 
30 years ago [21]. The algorithm uses standard linear and mixed integer linear programming, and 
it is based on elementary graph theory and important former results on the dense realizations of 
CRNs. The proposed method is also capable of determining if no dynamically equivalent weakly 
reversible structure exists for a given reaction network with a previously fixed complex set. 
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1 Introduction 



Chemical reaction networks (CRNs) are widely used for modeling chemical and biochemical pro- 
cesses in laboratory, industrial or natural environments [38]. Despite their simple mathematical 
structure, CRN models are capable of describing complex dynamical phenomena such as multiple 
steady states or oscillating behaviour that can be of key importance in understanding the oper- 
ation of the modeled technological or living mechanisms [27, 37, 31]. From the point of view of 
mathematical sciences and systems theory, CRN models and the related analysis methods are also 
interesting and well-utilizable in the study of the characteristic features of nonlinear dynamical 
systems [1, 6]. A special class of CRNs is the set of networks whose dynamics is governed by the 
mass action law (MAL). In this paper, we will deal with such systems. 

The foundations of chemical reaction network theory (CRNT) defining the basic notions and 
properties of reaction networks, were laid in the early papers and lecture notes [20, 14, 15, 16]. 
Since then, the development in the field has been continuous with significant contributions such as 
in [21, 13, 30, 18, 8, 7, 9, 29]. 

It has been known for long that CRNs with different structures and even with different set of 
chemical complexes may give rise to exactly the same set of differential equations describing the 
time-evolution of specie concentrations [20]. This fact is called macro- equivalence in [20] and it is 
sometimes referred to as the "fundamental dogma of chemical kinetics". Different CRNs producing 
the same ODEs will be called dynamically equivalent in this paper. Despite of the usefulness of 
dynamic equivalence in the analysis of CRNs, this property has only been studied to a rather 
limited extent in the literature. Conditions for the dynamical equivalence to a detailed balanced 
mechanism were first given in [2]. In [10], the authors gave necessary and sufficient conditions for 
the dynamic equivalence (also called confoundability) of CRNs that was extended with a special case 
in [32]. In [22], the notion of dynamical equivalence was significantly extended by the introduction 
of conjugate networks that mean CRNs with qualitatively the same dynamics where there is a 
well-defined mapping that takes trajectories of one system to into trajectories of the other. Also in 
[22], conditions for conjugacy are given when the above mentioned mapping is a nontrivial linear 
transformation. 

In [33] the terms dense and sparse realizations were introduced for dynamically equivalent 
reaction networks containing the maximal or minimal number of reactions with a fixed set of com- 
plexes. Additionally, a computational method based on mixed integer linear programming (MILP) 
was proposed to compute such realizations. In [36], important properties of dense realizations were 
described: it was shown that the structure (i.e. the unweighted directed graph) of a dense realiza- 
tion of a given CRN is unique, and the unweighted directed reaction graph of any other realization is 
the subgraph of the dense realization if the set of complexes is fixed. This is a key result that forms 
the basis of the algorithm presented in this paper. Furthermore, additional MILP based methods 
were also given in [36] for computing reversible CRN realizations, and realizations with the minimal 
number of complexes. It was shown in [35] that the computation of dynamically equivalent detailed 
balanced and complex balanced CRNs can be traced back to simple linear programming (LP). 

Many strong results of Chemical Reaction Network Theory apply to weakly reversible networks. 
The most well-known examples are probably the Deficiency Zero and Deficiency One Theorems 
relating complex composition, network structure and the qualitative properties of CRN dynamics 
[16]. Furthermore, it is known that for any weakly reversible network structure, there always exists 
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such a parametrization of the reaction rate coefficients that gives a complex balanced system [19]. 
Finally, the famous persistency conjecture, says that the dynamics of weakly reversible networks is 
persistent in the sense that no trajectory that starts in the positive orthant has an w-limit point 
on the boundary of . 

Similarly to numerous other key CRN properties such as deficiency, full reversibility, complex 
or detailed balance, weak reversibility is also not the property of the differential equations of the 
CRN models, but it is a realization property. Among the concluding open questions of [21] (where 
the authors gave necessary and sufficient conditions for polynomial ODEs to be kinetic together 
with a constructive proof) we can read the following: "VFe may look for a mechanism in a class of 
mechanisms with a given - chemically relevant - property. Such a property may be conservativity, 
(weak) reversibility, zero deficiency or just structural stability as well." Additionally, in [22] the 
following is written: "T/ie development of algorithms and computer software which can efficiently 
check for viable weakly reversible target networks . . . is a primary interest". Therefore it is of 
significant importance to work out a method that is able to construct a dynamically equivalent 
weakly reversible realization for a given CRN if it exists. 

Optimization tools based on e.g. linear /nonlinear programming or integer programming have 
been widely and successfully used in solving scientific and engineering problems [11, 12, 17]. Op- 
timization methods can be very useful tools for deciding the feasibility of complex problems and 
give feasible solutions (if they exist), even if the original problem is hard to treat analytically [5]. 
We will follow this approach in this paper. 

The structure of the paper is the following. In section 2 the applied models and relevant 
properties of CRNs will be summarized. Section 3 contains the main contribution, i.e. the algorithm 
for determining weakly reversible realizations of CRNs. Illustrative examples are shown in section 
4, while the conclusions can be found in section 5. 

2 Models and properties of chemical reaction networks 

2.1 Basic notions for the description of chemical reaction networks 

Following [15] and several other works, we will characterize CRNs with the following three sets. 
1.5 = {^1, . . . , Xn} is the set of species or chemical substances. 

2. C = {Ci, . . . , Cm} is the set of complexes. Formally, the complexes are represented as linear 
combinations of the species, i.e. 



where aij are nonnegative integers and are called the stoichiometric coefficients. 

3. TZ = {{Ci, Cj) I Ci,Cj € C, and Ci is transformed to Cj in the CRN} is the set of reactions. 
The relation {Ci,Cj) S TZ will be denoted as Cj — )> Cj. Moreover, a nonnegative weight, the 
reaction rate coefficient denoted by kij is assigned to each reaction Ci — >■ Cj. 



n 




(1) 



3 



The above characterization naturahy gives rise to the fohowing graph structure (often called 
'Feinberg-Horn- Jackson graph' or simply reaction graph) of a reaction network [16]. The weighted 
directed graph G = {V, E) of a reaction network consists of a finite nonempty set V of vertices and 
a finite set E of ordered pairs of distinct vertices called directed edges. The vertices correspond 
to the complexes, i.e. V = {Ci,C2, ■ ■ ■ Cm}, while the directed edges represent the reactions, i.e. 
(Cj ,Cj) G if complex Ci is transformed to Cj in the reaction network. The reaction rate coeffi- 
cients kij are assigned as positive weights to the corresponding directed edges Ci — > Cj in the graph. 
A set of complexes {Ci, . . . , Cfc} is called a linkage class of a CRN, if the complexes of the set are 
linked to each other in the reaction graph but not to any other complex. It is remarked that loops 
(i.e. directed edges that start and end at the same vertex) are not allowed in reaction graphs. A 
reaction network is called reversible, if whenever it contains the reaction Ci Cj , then the reverse 
reaction Cj <— Ci is also present in the CRN. A reaction network is called weakly reversible, if each 
complex in the reaction graph lies on at least one directed cycle (i.e. if complex Cj is reachable from 
complex Ci on a directed path in the reaction graph, then Cj is reachable from Cj on a directed 
path). 

A directed graph is called strongly connected if there exists a path from each vertex of the graph 
to any other vertex. A strongly connected component or simply strong component of a directed graph 
is a set of vertices such that the directed edges between them provide a directed path from any 
vertex of the set to any other vertex, and to which no additional vertex can be added (i.e. a 
maximal strongly connected subgraph) . Since a vertex is naturally reachable from itself through an 
empty path, strong components containing only one vertex will be called trivial strong components. 
Clearly, the vertices of the individual linkage classes of a weakly reversible CRN form the strong 
components of the reaction graph. 

Assuming mass- action kinetics, the following dynamical description will be used to describe the 
time-evolution of specie concentrations [15, 16]: 

± = Y-Ak- i^ix) (2) 

where Xi denotes the concentration of specie Xi . The ith column of Y contains the composition of 
complex Ci, i.e. [^]ij = ctji. The structure and parameters of the reaction graph are stored in the 
column conservation matrix Ak (also called the Kirchhoff matrix of the CRN) as follows 

t'^^^^^ -i if 

Finally, ip : M" i— ?• is a monomial-type vector mapping defined by 

n 

^,(x)=n^r", i=i,---,m (4) 
1=1 

2.2 Dynamically equivalent realizations 

As it has been mentioned in the introduction, reaction networks with different structures and/or 
parametrizations can give rise to the same kinetic differential equations. Therefore, we will call two 
reaction networks given by the matrix pairs {Y'^^\ A^^'^) and {Y^'^\ A^^^) dynamically equivalent, if 

y(l)^«^{l)(;r) = y(2)^{^2)^(2)^^^ ^ ^^^^^ ^ 



4 



where for i = 1,2, y(*) G R^x"^' have nonnegative integer entries. A, are vahd KirchhofF matrices, 
and 

mi 

V^f (x) = J]xf"^'='' ^ = 1,2, i = l,...,m. (6) 

k=l 

In this case, 

(y«^W) for i = 1,2 are cahed realizations of a kinetic vector field / (see, e.g. [21, 35] 

for more details). It is also appropriate to call {Y^^\ A^^}^) a realization of (Y^'^\ A^^^) and vica 
versa. 

3 The algorithm for computing weakly reversible realizations 

We recall that a dense realization of a CRN contains the maximal number of reactions (i.e. maximal 
number of nonzero off-diagonal elements in Af^) with a given stoichiometric matrix Y [33]. Based 
on the fact that with a given Y, all possible reactions are contained in the structurally unique dense 
realization [36], a straightforward idea is to try to find a dynamically equivalent weakly reversible 
mechanism starting from this superstructure. 

3.1 Basic principle of the algorithm 

Very shortly, the underlying principle of the presented algorithm is that it only removes (if possible) 
from the dense realization 

(i) edges that cannot be parts of any weakly reversible realization, 

(ii) edges the removal of which is necessarily implied by the deletion of edges belonging to set (i). 

The correct operation of the algorithm is based on the following two results. 

Rl If each strongly connected component of a directed graph G is contracted to a single vertex, 
the resulting directed graph is a directed acyclic graph [3] . (A directed graph is called acyclic 
if it has no nontrivial strongly connected subgraphs.) 

R2 The structure of the dense realization of any CRN is unique, and the directed unweighted 
graph of any CRN realization is a subgraph of the directed unweighted graph of the dense 
realization, if the set of complexes is fixed [36]. 

From R2 it follows that for obtaining a CRN superstructure including all possible structures, a 
dense realization must be computed. Directed edges between different strong components must be 
removed because they cannot lie on a directed cycle in any realization (Rl). If this is not possible, 
then there is no weakly reversible realization of the CRN. However, if the deletion is possible, it 
may imply the removal of additional reactions because of the linear kinetic constraints (see eqs. 
(9)-(12)). In general, this may impair the weak reversibility of the obtained network. In such 
a case, a new dense realization must be computed excluding the unnecessary edges identified in 
the previous step, and the procedure must be repeated until either a weakly reversible realization 
is found, or the deletion of undesired edges is no longer possible. In the latter case, no weakly 
reversible realization of the initial CRN exists with the given stoichiometric matrix Y. 
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3.2 Definition of input data structure and the necessary additional procedures 

We will assume that an initial CRN realization is given with the matrices {Y^^\ A^^^), and M = 
y(o) ^ T^liQ constraint set containing directed edges to be eliminated from the current realization 
is denoted as 

J<^ = {{Pi,qi),---,{Ps,qs)}, s<r. (7) 

where pi and qi denote the indices of the initial and terminal vertices of the ith edge, respectively, 
and r is the number of reactions in the CRN. 

Now, the following simple procedure will be defined for later use. 

L = FindCrossComponentEdges(^*fc") (8) 

The input of the procedure is the Kirchhoff matrix of a CRN, and the output is a set L con- 
taining the directed edges linking different strong components of the reaction graph. The strongly 
connected components of a directed graph can be determined in linear time using e.g. Kosaraju's, 
Tarjan's or Gabow's algorithm [3, 25]. Moreover, the examined CRN is weakly reversible if and 
only if it contains at least 2 reactions, and the output L of FindCrossComponentEdges is empty. 

In the following part of this subsection, appropriate modifications of the algorithm described in 
[33] are given, adapted to the current problem. 

3.2.1 Constraints corresponding to mass action dynamics 

The main characteristics of mass action dynamics are taken into consideration in the form of 
the following constraints (see also [33] for more details) in both of the procedures presented in 
subsections 3.2.2 and 3.2.3. 



Y-Ak = M (9) 

m 

Y,[Akh = 0, j = l,...,m (10) 

1=1 

[Ak\ij>0, i,j = l,...,m, i^j (11) 

[Ak]u<0, i = l,...,m (12) 



where Y and M are known, and the decision variables are the off-diagonal elements of Af^. It's 
easy to see that constraints (10)-(12) represent the fact that we are searching for a valid Kirchhoff 
matrix. 

3.2.2 Checking whether a set of reactions is removable from a CRN realization 

We call a set of reactions JC removable from a CRN realization, if there exists a dynamically 
equivalent CRN realization that does not contain the directed edges in /C. To check this, it is 
worth separately defining the following LP-based and thus polynomial time procedure to avoid 
unnecessary MILP computations that are known to be NP-hard. 
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The constraints (9)-(12) are completed with the following ones 

[■^klq^pi =0 for i = l,...,s (13) 

where {pi,qi) € /C as it is written in eq. (7). The feasibility of the linear constraints (9)-(12) and 
(13) can be checked by adding e.g. the following linear objective function to be minimized: 

m 

GLp{Ak)= Yl [^'^l^J (1^) 

hj = 1 

Clearly, eqs. (9)-(14) form a standard linear programming problem, the feasibility of which can be 
checked in polynomial time [11]. 

Based on the above, we will define the procedure to check whether a set /C of directed edges is 
removable from a CRN realization or not in the following way: 

Font = IsRemovable(y(*\ /C), (15) 

where the input data are as follows: {Y^^,Af) is a CRN realization, and /C is a constraint set 
of the form (7) containing the tail and head index pairs of the directed edges to be removed from 
(yi^^Af). The output Font is a boolean variable: its value is true if there exists a dynamically 

equivalent realization to (y^*\ A^*^) not containing the edges listed in K, and it is false if there 
does not exist such realization. 

3.2.3 Computing the dense realization excluding given directed edges 

For the solution of this subtask, the well-known results on the combination of propositional logic 
with mixed integer programming are applied [26, 4]. According to these, a propositional logic 
problem, where a statement must be proved to be true if a set of compound statements are also 
given, can be solved through a linear integer problem. 

If the procedure IsRemovable returns a true value, the dense realization of the CRN can be 
computed subject to the constraint that the edges listed in /C are excluded from it. To define the 
corresponding MILP problem, first we add exactly the same linear constraints contained in eqs. 
(9)-(13) as in the previous case. To make the forthcoming problem computationally tractable, we 
also introduce the following bounds for the decision variables 

[Ak]ij<Uij, Uij>0, i,j = l,...,m, i^j, (j,z)^/C (16) 
[Ak]ii>li, li<0, i = l,...,m (17) 

Here we are searching for such Ak that contains the maximal number of nonzero off-diagonal 
elements. For this, logical variables denoted by 5 are introduced and the following compound 
statements are constructed 

6ij = 1 ^ [Ak]ij > e, i,j = l,...,m, i / j, (j,i)^/C (18) 



7 



where the symbol "o" denotes "if and only if, and < e ^ 1 (i.e. elements of below e are 
treated as zero). Considering also (16), statement (18) can be translated to the following linear 
inequalities 

< [Aklij - e6ij , i,j = l,...,m, i ^ j, (j,z)^/C (19) 
< -[Ak]ij + Uijdij, i,j = l,...,m, i^j {j,i)^lC (20) 

Now it is possible to compute the realization containing the maximal number of reactions by 
maximizing the objective function 

m 

Gmilp{S)= (21) 

i 3, U, i) 

Finally, the following procedure is defined based on the MILP problem given by eqs. (9)-(13), 
(16)-(17), and (19)-(21). 

Aj. = FindConstrDenseRealization(y(*),^^*\/C), (22) 

where the input data set is the same as in the case of the procedure IsRemovable (see eq. (15) 
and the corresponding description). The output Aj, is a dense realization that does not contain 
the directed edges listed in the set /C. If the set /C is empty, then the procedure is the same that 
computes dense realizations and that was published in [33] . 

3.3 Formal description of the algorithm 

Now we are in the position to give the formal description of the procedure for determining weakly re- 
versible CRN realizations. The input data of the procedure is an initial CRN realization {Y^^\A^I!'^). 
The output is an m x m matrix that is the Kirchhoff matrix of the weakly reversible realization if 
there exists such, or a zero matrix if the procedure found no weakly reversible realizations. In the 
algorithm pseudocode, the auxiliary variable ExitCondition is a boolean storing the exit condition 
from the main loop. The complete pseudocode of the procedure called FindWeaklyReversibleRe- 
alization with common notations and keywords can be found in Table 1. 

3.4 Main properties of the algorithm 

The above described algorithm always finds a dynamically equivalent weakly reversible realization, 
if it exists. This clearly follows from the results Rl, R2 and the basic principles of the algorithm 
described in subsection 3.1. From these facts it also follows that the algorithm finds the 'densest' 
weakly reversible realization that structurally contains any other weakly reversible realizations 
(for illustration, see the results in subsection 4.1). An apparent drawback of the algorithm that it 
contains a step where a MILP problem has to be solved that is known to be NP-hard. However, the 
parallel (columnwise) implementation of the procedure FindConstrDenseRealization is possible 
as it is explained in [33] , that significantly extends the number of complexes that can be handled in 
practice. We remark that the immediate deletion of the columns/rows corresponding to the isolated 
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A°"*=FindWeaklyReversibleRealization(y(°),Af^) 

1 *:=0 € M™^™; ExitCondition:=fa.lse; 

2 Y := y(0). - ^(0). i^^^^-true; /C := {}; L := {}; 

3 while {ExitCondition= false) do 

4 begin 

5 if (/C 7^ {}) then Fotjj:=IsRemovable(y,^fc,/C); 

6 if (-Fout =true) then 

7 begin 

8 ylfc:=FindConstrDenseRealization(y,Afc,/C); 

9 L:=FindCrossComponentEdges(j4fc); 

10 if (L = {}) then ExitCondition:=true; Al^^:=Ak; 

11 else /C := /C U L; 

12 end 

13 else ExitCondition:=tzue; 

14 end 

15 return 



Table 1: Pseudocode of the algorithm for finding weakly reversible realizations 



complexes from matrices Y and after calling the procedure FindConstrDenseRealization 
is also possible, but this requires the renumbering of complexes. This is a technical detail of 
implementation and does not affect the principle or final output of the algorithm. (However, 
decreasing the number of optimization variables and constraints in each step in such a way might 
be required especially in the case of larger networks.) 

4 Examples 

The following examples were implemented in the MATLAB R13 computation environment using 
the YALMIP and the Multi-Parametric Toolboxes [24, 23]. The examples were run on a desktop 
PC with dual Intel Xeon 1.8GHz CPU and 2 Gigabytes of RAM. The implementation of dense 
realization computation was not parallel, therefore all the decision variables and constraints were 
put into a single optimization problem in procedure FindConstrDenseRealization. The strong 
components of the reaction graphs were identified using Kosaraju's algorithm [28]. 

4.1 Weakly reversible realizations of a simple irreversible network 

The simple network that can be seen in Fig. 1 a) was taken from [22] (Example 3). In [22], it 
is shown that for any positive e, the network has a possible weakly reversible realization with the 
structure shown in Fig. lb). (The computation and analysis of the parameters of the CRN in Fig. 
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Figure 1: a) Simple irreversible network from [22], b) Structure of one of its possible weakly 
reversible realizations determined in [22] 



lb can be found in [22].) The stoichiometric matrix of the network is 



Y 



The nonzero elements of the Kirchhoff matrix Ak € 



k\2,l 



1.5, \A 



fcj4,3 
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with e = 1.5 are 



1, [Akhfi 



1 . 



(23) 



(24) 



For this network, the algorithm described in section 3.3 and Table 1 works as follows. After the 
initialization steps, the dense realization containing all possible reactions with an empty constraint 
set /C is computed (line 8 of the pseudocode) . The Kirchhoff matrix of the dense realization is given 
by 



A 
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(25) 



The structure of the dense realization is shown in Fig. 2. The following steps can be followed using 
Fig. 3. The dense realization is not weakly reversible because there are edges between different 
strong components (line 9). The complexes of the single nontrivial strong component are indicated 
by boldface labels in Fig. 3. It is clear from the figure that the edges adjacent to the complexes 
Xi, 3X2, 3X1 + X2 are to be removed. These edges are drawn with dotted arrows in the figure. 
Thus, the constraint list is 

/C = {(1, 2), (1, 4), (1, 7), (3, 2), (3, 4), (3, 7), (5, 2), (5, 4), (5, 7), (6, 2), (6, 4), (6, 7)}. 

The next iteration of the algorithm (line 5) gives that these edges can be removed from the realiza- 
tion. Now a new dense realization is computed excluding edges in IC (line 8). The reactions of this 
dense realization are indicated by thick arrows in Fig. 3. Note that the constraint of excluding the 
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Figure 2: Structure of the dense realization of the reaction network shown in Fig. 1 a) 



reactions in /C from the network resulted in the removal of directed edges Xi + X2 — > Xi + 2X2, 
Xi + 2X2 2X1 + X2, Xi + 3X2 2X1 + X2, and Xi + X2 ^ Xi + 3X2, too, that were within 
the nontrivial strong component of the previous step. In this case, the resulting network remained 
weakly reversible (lines 9-10), so the algorithm can stop with success and return the determined 
dynamically equivalent weakly reversible CRN. The structure of the resulting network is shown in 
Fig. 4. The Kirchhoff matrix of the obtained realization is the following 
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The running time of the algorithm was lis using the hardware/software environment described in 
the beginning of section 4. 

It is interesting to note that the obtained weakly reversible realization is not complex balanced. 
(For convenience, we briefly define the notion of complex balance: A CRN realization {Y,Ai^) is 
called complex balanced if there exists a positive vector x* G for which ^^^(x*) = [20, 30, 
7].) However, using the polynomial-time algorithm described in [35], we can compute a complex 
balanced realization with 6 reactions in 0.1s, that is shown in Fig. 5. It is easy to see that the 
unweighted directed graphs of Figs. 1 b) and Fig. 5 are the proper subgraphs of the structure 
visible in Fig. 4. 
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Figure 3: Illustration of the operation of the algorithm on the example given in section 4.1 




Figure 4: The structure of the obtained dynamically equivalent weakly reversible CRN 




Figure 5: Complex balanced realization of the CRN described in section 4.1 
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4.2 Weakly reversible realization of a kinetic polynomial system 



The equations and initial CRN in this example were also used in [35] and [34] for illustrating other 
computation methods. We start from the polynomial ordinary differential equations given by 

Xi = — X1X2 + X3X4 — 2x1X2X3 
±2 = x\— X1X2 + 2x3X4 — 4x1X2X3 

X3 = —2x3 + X1X2 — xix|x3 + 2x4 (27) 
X4 = X1X2 — X3X4 + 4x1X2X3 — 3x| 

Using the inverse kinetic algorithm first published in [21] (see also [35] for a summary about 
the kinetic realizability of polynomial vector fields), a CRN shown in Fig. 6 can be constructed 
algorithmically that gives the dynamics (27). The numbering of complexes in the figure is the 
following: 

1 : 2X3, 2 : X3 + X4, 3 : Xi + 2X3, 4 : X2 + 2X3, 
5 : X3, 6 : Xi + X3 + X4, 7 : X2 + X3 + X4, 
8 : Xi + X2, 9 : Xi + 2X2 + X3, 10 : Xi, 11 : X2, 
l2:Xi + X2 + X4, 13 : Xi + X2 + X3, 14 : 2X2 + X3, 
15 : Xi + 2X2, 16 : Xi + 2X2 + X3 + X4, 

17 : 3X4, 18 : X3 + 3X4, 19 : 2X4 (28) 

The dense realization of the network contains 80 reactions, therefore it is not shown in a figure. 
For the sake of completeness, the list of the reactions (i.e. weighted directed edges) of the dense 
realization in the form [source vertex number, destination vertex number, rate coefficient) is given 
below: 

(5, 1, 0.5), (7, 1, 0.1), (11, 1, 0.1), (15, 1, 0.8), (1, 2, 0.1), (3, 2, 0.1), (5, 2, 0.1), (7, 2, 0.1), (12, 
2, 0.3), (1, 3, 0.1), (5, 3, 0.1), (7, 3, 0.1), (12, 3, 0.1), (1, 4, 0.1), (3, 4, 0.4), (5, 4, 0.1), (7, 4, 0.1), 
(11, 4, 0.1), (3, 5, 0.1), (7, 5, 0.1), (11, 5, 0.1), (15, 5, 0.1), (3, 6, 0.1), (5, 6, 0.1), (7, 6, 0.1), (1, 7, 
0.1), (3, 7, 0.1), (5, 7, 0.2), (12, 7, 0.3), (1, 8, 0.1), (3, 8, 0.1), (5, 8, 0.3), (7, 8, 0.3), (11, 8, 1.2), 
(1, 9, 0.1), (5, 9, 0.1), (7, 9, 0.1), (11, 9, 0.2), (1, 10, 0.5), (3, 10, 0.8), (5, 10, 0.1), (7, 10, 0.1), (12, 
10, 0.1), (3, 11, 0.1), (5, 11, 0.1), (7, 11, 0.1), (1, 12, 0.1), (3, 12, 0.1), (5, 12, 0.1), (7, 12, 0.1), (1, 
13, 0.1), (3, 13, 0.1), (5, 13, 0.1), (7, 13, 0.1), (11, 13, 0.1), (15, 13, 0.1), (1, 14, 0.1), (3, 14, 0.1), 
(5, 14, 0.1), (7, 14, 0.1), (12, 14, 0.1), (3, 15, 0.1), (5, 15, 0.1), (7, 15, 0.7), (11, 15, 0.1), (5, 16, 
0.25), (7, 16, 0.3), (11, 16, 0.7), (15, 16, 0.2), (3, 17, 0.1), (5, 17, 0.1), (7, 17, 0.1), (3, 18, 0.1), (5, 
18, 0.1), (7, 18, 0.4), (3, 19, 0.1), (5, 19, 0.1), (7, 19, 0.1), (11, 19, 0.1), (15, 19, 0.1). 

The dense realization has 13 strong components. Out of these, there is only one nontrivial 
strongly connected component containing the vertices 1, 3, 5, 7, 11, 12, and 15. After identifying 
all directed edges linking different strong components, we obtain that the following 53 edges given 
in the form [source vertex number, destination vertex number) should be deleted from the dense 
realization in the next step of the algorithm: 

(1,2), (3,2), (5,2), (7,2), (12,2), (1,4), (3,4), (5,4), (7,4), (11,4), (3,6), (5,6), (7,6), (1,8), (3,8), (5,8), 
(7,8), (11,8), (1,9), (5,9), (7,9), (11,9), (1,10), (3,10), (5,10), (7,10), (12,10), (1,13), (3,13), (5,13), 
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Figure 6: CRN realizing the dynamics of eq. (27). Only reaction rates different from 1 are indicated. 



(7,13), (11,13), (15,13), (1,14), (3,14), (5,14), (7,14), (12,14), (5,16), (7,16), (11,16), (15,16), (3,17), 
(5,17), (7,17), (3,18), (5,18), (7,18), (3,19), (5,19), (7,19), (11,19), (15,19). 

All the above listed edges were possible to remove from the dense realization (their removal 
implied the deletion of 12 additonal reactions), and the resulting constrained dense realization is 
shown in Fig. 7. It is clear from the figure that edges adjacent to vertices no. 11 and 12 have to be 
removed in the following step. This final step is illustrated in Fig. 8, where the meaning of the line 
types is the same as in the case of Fig. 3. With the removal of edges (5,11), (5,12), (7,11), (7,12), 
the directed edges (5,1), (7,1), (7,3), (5,3), (7,5) and (5,15) were also deleted, and the resulting 
weakly reversible realization (of deficiency 0) is shown in Fig. 9. The total running time of the 
algorithm was 80.5s. 




Figure 7: CRN realization of eq. (27) containing the vertices of the only nontrivial strong component 
of the dense realization 
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5 Conclusions 



A numerical algorithm for finding dynamically equivalent weakly reversible realizations of chemical 
reaction networks was proposed in this paper. The algorithm computes a maximal superstructure 
that contains all other possible weakly reversible structures as proper subgraphs, and it is able 
to determine if no weakly reversible realization exists. The method uses linear and mixed integer 
linear programming steps and it is based primarily on the computation and properties of dense 
realizations published in [33, 36]. The computationally critical MILP step of the algorithm can be 
implemented parallelly. The operation of the algorithm has been illustrated on numerical examples. 
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